Observation of shock waves in a large Bose-Einstein condensate 
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We observe the formation of shock waves in a Bose-Einstein condensate containing a large number 
of sodium atoms. The shock wave is initiated with a repulsive, blue-detuned light barrier, intersecting 
the BEC, after which two shock fronts appear. We observe breaking of these waves when the size 
of these waves approaches the healing length of the condensate. At this time, the wave front splits 
into two parts and clear fringes appear. The experiment is modeled using an effective ID Gross- 
Pitaevskii-like equation and gives excellent quantitative agreement with the experiment, even though 
matter waves with wavelengths two orders of magnitude smaller than the healing length are present. 
In these experiments, no significant heating or particle loss is observed. 



I. INTRODUCTION 

The realization of Bose-Einstein condensates (BECs) 
in dilute atomic gases [l[ provides the opportunity for 
the study of non-linear matter wave dynamics. Many 
experiments on both the statics and dynamics of BECs 
have shown that experiments can often be modeled ac- 
curately by solving the mean-field Gross-Pitaevskii equa- 
tion (GPE) @, i, 1, |, i, 0, I, B M HI E3, 0- For 
example it has successfully been used to model experi- 
ments on interferometry soliton formation El, El, 
four wave mixing atom laser outcoupling IToLTll j . 
sound propagation and superfluidity [2|] to name a few. 
However, the Gross-Pitaevskii equation is not an exact 
description of a Bose-Einstein condensate, but instead it 
is expected to be a good approximation for condensates 
that contain a relatively large number ofparticles and 
are not undergoing violent dynamics [H], [l5[ . Thus it is 
of interest to experimentally probe the limits of validity 
of the GPE in describing the dynamics of Bose-Einstein 
condensates. 

The well-known 'Bosenova' experiment of Donley et 
al. [13] is one situation where the validity of the GPE 
might be questioned, combining a small number of atoms 
with violent dynamics. The experiment began with a 
near ideal 85 Rb condensate of a few thousand atoms in 
its ground state before the atomic scattering length was 
manipulated using a Feshbach resonance to being attrac- 
tive. The BEC was observed to collapse, emitting high 
energy jets and bursts of atoms. A number of compu- 
tational studies have shown good qualitative agreement 
with the experimental observations [18l [l9[ |20| . However, 
careful quantitative studies have indicated that at the 
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most basic level there is quantitative disagreement with 
the experimental data. Savage et al. showed that the 
experimentally measured collapse time is typically 25% 
shorter than predicted by simulations of the GPE [2lj . 
Going beyond mean-field theory and including the first 
order quantum corrections to the dynamics was shown to 
make little difference to the numerical results [22|, |23[ . A 
second group of experiments by the same group found ev- 
idence for the formation of repulsive bright solitary waves 
in the condensate collapse |24 |. However, modelling this 
experiment with the GPE is unable to reproduce this ex- 
perimental finding [251 ]. 

Another example where the GPE has successfully mod- 
eled results of violent experiments is in the formation 
of quasi- ID bright solitons as reported in Refs. [2(1 127|]. 
Both of these experiments involve a sudden change to 
an attractive scattering length. It has been shown Q 
that the GPE can be used to successfully model these 
experiments. However, in this case, the three-body re- 
combination, which is not included in the GPE deriva- 
tion, plays a nontrivial role in the dynamics, and the GPE 
must be modified by the addition of phenomenological 
damping terms in order to agree with these experiments. 

Another situation exhibiting violent dynamics in the 
solution of the GPE without dissipation is in the gener- 
ation of shock waves. Damski [16] has calculated the ID 
GPE dynamics following the introduction and sudden re- 
moval of an attractive dimple potential in the centre of a 
harmonically trapped elongated BEC. The localised den- 
sity bulge splits into two pulses, which propagate towards 
the ends of the condensate. Due to the density depen- 
dent speed of sound in the system, the center of the pulse 
catches up to the front, creating a shock, and at this point 
the calculations develop a strong fringe pattern with a 
spacing of the order of the healing length corresponding 
to classical wave breaking. However, Damski speculated 
that in a physical system the Gross-Pitaevskii equation 
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would become invalid at this point and this would not be 
observed in an experiment. 

Recently, a number of experiments have been per- 
formed that have observed phenomena related to wave 
breaking in a Bose-Einstein condensate. The first exper- 
iments were by Dutton et al. @ . They used the slow light 
mechanism to create a defect in the condensate which was 
much narrower than the healing length. This defect cre- 
ated dark solitons which shed high frequencies traveling 
at different velocities. The wave front of the propagat- 
ing solitons eventually became curled and decayed into 
vortex pairs. 

Simula et al. [3l| blasted a hole in a rapidly rotating 
oblate BEC using a repulsive dipole potential and found 
qualitative agreement between their experimental obser- 
vations and numerical calculations of the GPE. This was 
followed by Hoefer et al. who performed similar exper- 
iments on a stationary and slowly rotating oblate BEC, 
and made the connection to dispersive shock-waves in 
the GPE [32]. They found that this wave breaking phe- 
nomenon can be described by the Gross-Pitaevskii equa- 
tion [33]. However, to obtain quantitative agreement be- 
tween the simulations and the experimental observations 
they found that they had to use the laser width as a fit- 
ting parameter. In one situation agreement was obtained 
with a value 1.5 times larger than in the experiment, but 
in another it was required to be half as large, suggesting 
that this is not a systematic error but something more 
fundamental. 

The purpose of this paper is to generate shock waves in 
a large number, elongated BEC and to perform a quan- 
titative comparison with simulations of an effective ID 
GPE as a test of its validity in extreme conditions. We 
observe no indication of heating or particle loss, and con- 
clude that even at these relatively high energy scales the 
GPE is valid for the time of the experiment. Soon af- 
ter the experimental data for this work was taken, re- 
lated work was reported by Chang et al. (34| . They also 
found good agreement with the GPE. However, these ex- 
periments were performed at significantly lower energy 
scales, where agreement with the GPE is expected, and 
hence they did not perform a detailed investigation of the 
fringe spacing and possibility of heating. 



II. EXPERIMENT 

In our experiment we begin with an almost pure 
sodium BEC containing up to 100 million atoms in the 
|F, mp) = |1,— 1) hyperfine state. The atoms are mag- 
netically confined in a clover leaf trap with harmonic 
trapping frequencies of 97 Hz in the radial direction and 
3.9 Hz in the longitudinal direction [36[. For the conden- 
sate, this corresponds to a chemical potential of 214 nK 
and Thomas-Fermi widths of 20 /am and 506 /am in the 
radial and axial directions, respectively. To create the 
initial disturbance we turn on a repulsive, blue-detuned 
laser beam focussed in the middle of the BEC, intersect- 



ing the cloud in the longitudinal direction. The focus 
has a waist of 90 ± 11 /im (1/e 2 ) and a wavelength Al, 
tunable from 567 to 584 nm, where the power can be 
switched within 10 fis. At t = the blue-detuned laser is 
turned on suddenly with the magnetic trap still on and 
after variable times the cloud is imaged in the trap. Note 
that the waist of the laser focus is much smaller than the 
axial width of the condensate and much larger than the 
radial widths, suggesting that most of the dynamics will 
occur along the axial dimension. 



A. Case 1 

For the first set of experiments we make use of phase- 
contrast in-trap imaging (3?J . The detuning of the imag- 
ing laser is chosen such that a 2tt phase shift is accumu- 
lated for the largest atom cloud density. The resolution is 
3 x 3 jam and the images are recorded with a CCD camera 
(Apogee AP1E). Figure [1] shows the in-trap CCD images 
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FIG. 1: Phase contrast images of the BEC in the trap for (A- 
D) t = 2, 5, 8 and 15 ms after turning on the blue-detuned 
laser beam. The image size is 1200 x 300 /xm. 

for a BEC with 50 million atoms, corresponding to a 
chemical potential \i = 162 nK, imaged after t = 2,6,8 
and 15 ms respectively. The power for the repulsive laser 
beam used is 78 mW with a wavelength of 579 nm, result- 
ing in a repulsive barrier of 12.6 /iK. For all experiments 
the scattering of photons from the laser can be neglected. 
The first image shows the splitting of the cloud induced 
by the switch on of the repulsive barrier. The density 
profile of the BEC results in a gradient in the speed of 
sound of the condensate v c as v c ~ y/n. As a conse- 
quence, as the wave fronts induced by the lasers travel 
from the centre towards the edge of the condensate width 
the back of the pulse catches up to the front, resulting in 
steepening of the wave (self-steepening). This is clearly 
observed in the second image. The last two images show 



that the wave breaks into two parts after the maximum 
density gradient has been reached. 



B. Case 2 

As the size of the fringes in the data set above was 
much smaller than the resolution of the imaging system, 
the experiment was repeated and imaged using a 70 ms 
time of flight expansion to enlarge the detailed features 
of the condensate. This free expansion time was chosen 
in such a way that the full CCD array is used to obtain 
the maximum resolution. Figure [2] shows the images for 
a smaller condensate of 18 million atoms for an evolution 
time in the trap of 0, 1, 2, 6, 13 and 25 ms respectively. In 
this situation we use a repulsive laser beam with a power 
of 21 mW and a wavelength of 569 nm, corresponding to 
a repulsive potential barrier with a height of 0.43 /iK. 

The images in Fig. [2] indicated that it takes longer for 
the build-up of the wave front for a smaller BEC and a 
weaker repulsive barrier. It can be seen that the shock 
front splits into two parts, where pronounced fringe pat- 
terns appear. Fourier analysis on the last image shows a 
fringe size of 6 ± 1/im. Integration over the entire pixel 
array shows a 2.7% fluctuation in the number of atoms 
over a evolution time of 30 ms, suggesting that losses 
due to any potential heating are negligible. Frames (E) 
and (F) show that there is a slight curvature of the wave 
fronts, suggesting that the outer edge of the wave front 
has travelled faster than the center. 



III. NUMERICAL SIMULATIONS 
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FIG. 2: Absorption images of an 70 ms expanded BEC for (A- 
F) t = 0, 1, 2, 6, 13 and 25 ms after turning on the repulsive 
barrier. The size of the images is 2304 x 1536 fim. 



We now turn to numerical modelling these experi- 
ments. Ideally we would do so with a fully three- 
dimensional simulation of the Gross-Pitaevskii equation. 
However, even making use of the cylindrical symmetry we 
find that the simulation parameters for such a large BEC 
are too demanding for numerical simulation with the 
computational resources we have available. To proceed 
we have made use of the fact that the system has cylin- 
drical symmetry and a high aspect ratio, and reduced 
the problem to simulating the one-dimensional nonpoly- 
nomial Gross-Pitaevskii equation as derived by Salasnich 
et al. [38[ (NPGPE). This assumes a variational gaussian 
profile for the wave function in the radial direction whose 
width, cr(z), is dependent on the local one dimensional 
density. This is to allow for the 'bulging' in the radial di- 
rection due to the mean field interaction to be included. 
The result is a ID Gross-Pitaevskii equation where the 
nonlinear interaction term is a function of the local den- 
sity. In the case of a harmonic trapping potential in the 
radial direction with trapping frequency oj±_ , the effective 



ID equation of motion is 
with 

a 2 (z) = -^- v /l + 2a s |^( 2 )| 2 . (2) 
mco± 

V(z) represents the external potential along the z direc- 
tion, Usd = 4:7rh 2 a s /m is the 3-d atom-atom interaction 
strength, and a s is the s-wave scattering length. We solve 
this equation numerically with no free parameters. 

A. Case 1 

We used Eq. (pQ) and Eq. ([2j) to simulate the experimen- 
tal conditions used in case 1. We used an initial condition 



of 50 million atoms in the ground state of a harmonic trap 
of trapping frequencies uo z = 2tt x 3.9 Hz and co± = 2tt x 97 
Hz. The ground state was found by imaginary time prop- 
agation. A repulsive potential due to the laser beam was 
then switched on. The laser power was 49 mW, and de- 
tuned by 23 nm from resonance. We assumed that the 
beam profile was gaussian with a width of 90 /im. This 
corresponds to a repulsive potential height of 0.87 jaK. 

Figure [3] shows a comparison of theoretical calculations 
with experimental data for these parameters. The left 
column (a-d) shows the density profile of the condensate, 
integrated along the radial direction, and the right col- 
umn (e-h) shows the theoretical calculations. The images 
are for (a) and (e) 11 ms, (b) and (f) 14 ms, (c) and (g)21 
ms (d) and (h) 32 ms after the repulsive potential was 
switched on. 
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FIG. 3: Comparison of theoretical calculations with experi- 
mental data. The left column ((a)-(d)) shows the density pro- 
file of the condensate, integrated along the radial direction, 
and the right column ((e) -(h)) shows the theoretical calcula- 
tions. The images are for (a) and (e): 11 ms, (b) and (f): 14 
ms, (c) and (g): 21 ms (d) and (h): 32 ms after the repul- 
sive potential was switched on. The laser power was 49 mW, 
which corresponds to a repulsive potential height of 0.87 /iK 



to the inverse of the healing length 
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where U(z) = Usd/27ra 2 (z) is the effective 1-D interac- 
tion strength. In (a) and (d), \g(z)\ <C k w b(z), and the 
density profile of the pulse remains smooth. In (b) and 
(e), \g(z)\ becomes comparable to k w b(z), and a shock 
front begins to form. In (c) and (f), \g(z)\ ^> k w b, and 
high frequency oscillations are visible in the density pro- 
file. 



t = 14 ms 



300 r 
250 
200 
150 
100 
50 



(a) 



t = 15 ms 
(b) 



16 ms 



'(c) 




1^ 



(d) 



-2 [ 

0.26 



(e) 



(f) 








] 









0.27 0.28 0.29 0.27 0.28 

Z (mm) i 



0.29 0.30 0.28 0.29 0.30 0.31 

(mm) Z (mm) 



FIG. 4: Top row: Close up of the density profile of the leading 
edge of the pulse for (a): t — 14 ms, (b): t — 15 ms, and (c): 
t — 16 ms, for the same parameters as [3] Bottom row: The 
relative density gradient, g(z) (solid line) and the inverse of 
the healing length ±k w b(z) (dashed line) for the cases in (a), 
(b), and (c) respectively. In (d), \g(z)\ <C k w b(z), so the 
pulse propagates smoothly. In (e), g(z) becomes comparable 
to k w b(z), and a shock front begins to form. In (f), \g(z)\ ^> 
k w b, and high frequency oscillations are visible in the density 
profile. 

We found similar agreement between theoretical re- 
sults and experimental data for repulsive barrier heights 
of 0.39 /iK and 3.27 /iK (corresponding to data presented 
in figure (pQ)), which we have not presented here. 



The experimental images show that we get two pulses 
propagating outwards. The leading edges of these pulses 
gets steeper due to the density-dependent speed of sound 
[161 ]. The theoretical plots show agreement with the ex- 
perimental data in terms of the position of the pulses 
and the steepening of the density. However, they show 
extra detail which the experimental image is insufficient 
to resolve. As the pulses steepen the relative density 
gradient become comparable to the healing length, and 
a shock wave forms. This causes an abrupt change in 
the density profile of the smooth pulse, resulting in rapid 
density oscillations in front and behind the pulse. Figure 
(jl]) compares the relative density gradient 

dn(z) 



B. Case 2 

In the second set of data we made use of time-of-flight 
imaging and were able to observe features that were not 
resolvable in trap. To make a comparison with this ex- 
perimental data, we have modeled the expansion of the 
condensate in the radial direction during time of flight 
combined with continuing the ID simulation of the dy- 
namics of the wave function in the axial direction. We are 
unable to use the NPGPE in this case, as it is unable to 
model the expansion in the radial direction. Instead, we 
obtained a 1-D set of equations by assuming a transverse 
radial width of the condensate which was independent of 
z, which allowed us to scale our nonlinearity to one di- 
mension using Uid = Usd/itR 2 . We chose this width pa- 
rameter R by comparing the evolution while still in the 
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trap to the evolution of the NPGPE with no free param- 
eters, and chose the transverse width parameter which 
gave the best agreement in the results. Figure ([5]) shows 
the comparison between the solution to the GPE and the 
NPGPE after 3 ms of evolution, during the trapped phase 
of the evolution. For the GPE simulation, we set R = 12 
/im, which gave the best agreement with the NPGPE. 
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FIG. 5: Comparison of the density profiles resulting from 
GPE (blue solid line) and NPGPE (dashed black line) sim- 
ulations. The condensate was left to evolve in trap for 3 ms 
with a 0.88 /iK repulsive barrier. 

During the expansion phase, we assumed that the 
transverse width expanded according to the analytic re- 
sult derived by Castin and Dum [39[ for the self-similar 
expansion of a Thomas-Fermi condensate. According 
to this result, the transve rse width evolves according to 
R(t) = R(0)^l + (u; ± t) 2 . As we have neglected to in- 
clude three dimensional effects in our simulation, we do 
not expect perfect agreement with the experimental data. 

Figure [6] shows the a density slice from the experimen- 
tal data from Figure [2] (integrated in the y direction for 
100 pixels) compared to theoretical simulations. There is 
excellent agreement between these results. In both cases, 
high contrast fringes are observed when the back edge of 
the pulse catches up with the front edge. We find that 
the wave breaking occurs after a much shorter time com- 
pared to the case I, because the expansion in the radial 
direction leads to a rapid increase in the healing length. 

There is considerable agreement between our simula- 
tion with no free parameters, and the experimental data, 
indicating that it is valid to simulate such a violent sys- 
tem with the GPE. The main cause of discrepancy be- 
tween the experimental images and the theoretical cal- 
culations is most likely the uncertainty in the size of the 
waist of the blue detuned beam (90 ± 11 /im). By adjust- 
ing the value of the waist slightly in the calculations, we 
found better agreement with the experimental images, 
with the positions and widths of the wave packets giv- 
ing best agreement for a waist of 92 /im. There is also 
a slight discrepancy in the spatial frequency of the inter- 
ference fringes for the experimental images and the sim- 
ulation, the cause of which is unknown. It is difficult to 
accurately determine the spatial frequency in the experi- 
mental images, due to the 3 /im pixel size causing spatial 
aliasing. The drop in observed fringe contrast in the ex- 
perimental images is partly due to this spatial aliasing 
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FIG. 6: Density slice (left column), and theoretical 1-D GPE 
simulation (right column) of the density profile of the conden- 
sate after expansion for 69 ms after evolving in the trap for 
various times. The condensate was held in the trap with the 
optical dipole on potential for (a) and (f): 0.5 ms, (b) and 
(g): 1.1 ms, (c) and (h): 1.5 ms, (d) and (i):2.0 ms, (e) and 
(j):3.0ms. Parameters: N = 1.8 X 10 7 , waist = 90 /mi, P = 21 
mW, detuning =10 nm, corresponding to a barrier height of 
0.88 fj, K. 

effect, but most likely is the result of a slight misalign- 
ment of the imaging axis with respect to the axis of the 
fringes in the 3-D experiment. The asymmetry is due to 
a slight misalignment of the dipole barrier with the cen- 
tre of the trap. Figure © shows a comparison between 
the experimental data, and the results from the simula- 
tion convolved with a gaussian of width 3 /im, to emulate 
the effect of the 3 /im resolution. The fringe contrast is 
reduced as a result of the convolution, especially towards 
the back end of the pulse where the fringe separation is 
smaller. However, this effect is not enough to describe 
the lack of fringe contrast in the experimental images at 
the front of the pulse. One might suspect that the re- 
duced fringe visibility is due to heating of the condensate 
atoms. However, no thermal fraction was observed. This 
was validated by repeating the experiment for different 
expansion times. Inspection of Figure [2] shows that the 
wave front of the leading edge of the pulse is slightly 
curved. As this curvature will also be present in the di- 
rection of imaging, this will also contribute to the loss of 
fringe visibility. Given these discrepancies, it is clear that 
the solutions of the GPE correctly capture the physical 
processes determining the evolution of the system. 

IV. CONCLUSION 

We have generated shock waves in a large number elon- 
gated Bose-Einstein condensate by suddenly splitting it 
with a blue-detuned optical dipole potential. We ob- 
serve no significant particle loss or heating over these 
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FIG. 7: Close up view of the density slice for the same con- 
ditions as Figure [6] (e) and (j). (a) shows the raw simulation 
data, (b) shows the same data convolved with a 3 \i m gaussian 
to emulate the finite pixel size of the camera, and (c) shows 
the experimental data. The discrepancy in the positions of 
the experimental and theoretical density distributions is most 
likely due to the uncertainty in the waist of the blue-detuned 
beam. 



time scales, and find excellent agreement with the predic- 



tions of our simulations of the ID non-polynomial Gross- 
Pitaevskii equation. We used a method to describe the 
expansion of the elongated condensate while simulating 
the continuing dynamics in the long direction, and again 
find excellent agreement with simulations. Our study 
provides further evidence that the Gross-Pitaevskii equa- 
tion can be an excellent approximation to the dynamics 
of condensates even in situations exhibiting violent dy- 
namics and low dissipation. 
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